## Function simulating IPM model (end of the code) and Bootstrapping code ## parms=read.csv("IPM.parameters.csv") Nrep=500 result=data.frame(species=c(rep("BT",6),rep("TT",6)),herb=rep(c(rep("ambient",3),rep("reduced",3)),2), comp=rep(c("low","medium","high"),4), lambda=0, LCI=0, UCI=0) #BT, ambient herbivory, low competition resBT1=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH p.vec[3] = parms$ResVar.nobolt[1] #Residual variance of regression line p.vec[6] = parms$ResVar.bolt[1] #Residual variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) p.vec[13] = parms$size.dist.resse[1] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[1],parms$growth.intercept.nobolt.se[1]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[1],parms$growth.slope.nobolt.se[1]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[1],parms$growth.intercept.bolt.se[1]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[1] ,parms$growth.slope.bolt.se[1]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[1],parms$surv.intercept.se[1]) p.vec[8] = rnorm(1,parms$surv.slope[1],parms$surv.slope.se[1]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[1],parms$bolt.intercept.se[1]) p.vec[10] = rnorm(1,parms$bolt.slope[1],parms$bolt.slope.se[1]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[1],parms$recruit.se[1]) #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) p.vec[12] = rnorm(1, parms$size.dist.mean[1],parms$size.dist.unconditional.se[1]) #Average seedling size #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[2],parms$zero.seed.intercept.se[2]) p.vec[15] = rnorm(1,parms$zero.seed.slope[2],parms$zero.seed.slope.se[2]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[2],parms$seed.intercept.se[2]) p.vec[17] = rnorm(1,parms$seed.slope[2],parms$seed.slope.se[2]) resBT1[n]=ParBoot(p.vec) } pick = (result$species =="BT") & (result$herb == "ambient")& (result$comp=="low") result[pick,]$lambda=median(resBT1) result[pick,]$LCI=quantile(resBT1,0.05) result[pick,]$UCI=quantile(resBT1,0.95) ########################################################## #BT, ambient herbivory, med competition resBT2=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[2] #Residual variance of regression line #GROWTH - BOLTING p.vec[6] = parms$ResVar.bolt[2] #Residual variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[2] #Average seedling size p.vec[13] = parms$size.dist.resse[2] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[2],parms$growth.intercept.nobolt.se[2]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[2],parms$growth.slope.nobolt.se[2]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[2],parms$growth.intercept.bolt.se[2]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[2] ,parms$growth.slope.bolt.se[2]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[2],parms$surv.intercept.se[2]) p.vec[8] = rnorm(1,parms$surv.slope[2],parms$surv.slope.se[2]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[2],parms$bolt.intercept.se[2]) p.vec[10] = rnorm(1,parms$bolt.slope[2],parms$bolt.slope.se[2]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[2],parms$recruit.se[2]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[2],parms$size.dist.unconditional.se[2]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[2],parms$zero.seed.intercept.se[2]) p.vec[15] = rnorm(1,parms$zero.seed.slope[2],parms$zero.seed.slope.se[2]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[2],parms$seed.intercept.se[2]) p.vec[17] = rnorm(1,parms$seed.slope[2],parms$seed.slope.se[2]) resBT2[n]=ParBoot(p.vec) } pick = (result$species =="BT") & (result$herb == "ambient")& (result$comp=="medium") result[pick,]$lambda=median(resBT2) result[pick,]$LCI=quantile(resBT2,0.05) result[pick,]$UCI=quantile(resBT2,0.95) ########################################################## #BT, ambient herbivory, high competition resBT3=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[3] #Residual variance of regression line #GROWTH - BOLTING p.vec[6] = parms$ResVar.bolt[3] #Residual variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[3] #Average seedling size p.vec[13] = parms$size.dist.resse[3] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[3],parms$growth.intercept.nobolt.se[3]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[3],parms$growth.slope.nobolt.se[3]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[3],parms$growth.intercept.bolt.se[3]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[3] ,parms$growth.slope.bolt.se[3]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[3],parms$surv.intercept.se[3]) p.vec[8] = rnorm(1,parms$surv.slope[3],parms$surv.slope.se[3]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[3],parms$bolt.intercept.se[3]) p.vec[10] = rnorm(1,parms$bolt.slope[3],parms$bolt.slope.se[3]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[3],parms$recruit.se[3]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[3],parms$size.dist.unconditional.se[3]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[3],parms$zero.seed.intercept.se[3]) p.vec[15] = rnorm(1,parms$zero.seed.slope[3],parms$zero.seed.slope.se[3]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[3],parms$seed.intercept.se[3]) p.vec[17] = rnorm(1,parms$seed.slope[3],parms$seed.slope.se[3]) resBT3[n]=ParBoot(p.vec) } pick = (result$species =="BT") & (result$herb == "ambient")& (result$comp=="high") result[pick,]$lambda=median(resBT3) result[pick,]$LCI=quantile(resBT3,0.05) result[pick,]$UCI=quantile(resBT3,0.95) ################################################################# #BT, reduced herbivory, low competition resBT4=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[4] #Residula variance of regression line #GROWTH - BOLTING p.vec[6] = parms$ResVar.bolt[4] #Residula variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[4] #Average seedling size p.vec[13] = parms$size.dist.resse[4] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[4],parms$growth.intercept.nobolt.se[4]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[4],parms$growth.slope.nobolt.se[4]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[4],parms$growth.intercept.bolt.se[4]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[4] ,parms$growth.slope.bolt.se[4]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[4],parms$surv.intercept.se[4]) p.vec[8] = rnorm(1,parms$surv.slope[4],parms$surv.slope.se[4]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[4],parms$bolt.intercept.se[4]) p.vec[10] = rnorm(1,parms$bolt.slope[4],parms$bolt.slope.se[4]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[4],parms$recruit.se[4]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[4],parms$size.dist.unconditional.se[4]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS #all plants carried seeds p.vec[14] = 999 p.vec[15] = 999 #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[5],parms$seed.intercept.se[5]) p.vec[17] = rnorm(1,parms$seed.slope[5],parms$seed.slope.se[5]) resBT4[n]=ParBoot(p.vec) } pick = (result$species =="BT") & (result$herb == "reduced")& (result$comp=="low") result[pick,]$lambda=median(resBT4) result[pick,]$LCI=quantile(resBT4,0.05) result[pick,]$UCI=quantile(resBT4,0.95) #################################################################### #BT reduced herbivory medium competition resBT5=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[5] #Residual variance of regression line #GROWTH - BOLTING p.vec[6] = parms$ResVar.bolt[5] #Residual variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[5] #Average seedling size p.vec[13] = parms$size.dist.resse[5] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[5],parms$growth.intercept.nobolt.se[5]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[5],parms$growth.slope.nobolt.se[5]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[5],parms$growth.intercept.bolt.se[5]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[5] ,parms$growth.slope.bolt.se[5]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[5],parms$surv.intercept.se[5]) p.vec[8] = rnorm(1,parms$surv.slope[5],parms$surv.slope.se[5]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[5],parms$bolt.intercept.se[5]) p.vec[10] = rnorm(1,parms$bolt.slope[5],parms$bolt.slope.se[5]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[5],parms$recruit.se[5]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[5],parms$size.dist.unconditional.se[5]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[2],parms$zero.seed.intercept.se[2]) p.vec[15] = rnorm(1,parms$zero.seed.slope[2],parms$zero.seed.slope.se[2]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[2],parms$seed.intercept.se[2]) p.vec[17] = rnorm(1,parms$seed.slope[2],parms$seed.slope.se[2]) resBT5[n]=ParBoot(p.vec) } pick = (result$species =="BT") & (result$herb == "reduced")& (result$comp=="medium") result[pick,]$lambda=median(resBT5) result[pick,]$LCI=quantile(resBT5,0.05) result[pick,]$UCI=quantile(resBT5,0.95) ############################################################## #BT reduced herbivory, high competition resBT6=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[6] #Residula variance of regression line #GROWTH - BOLTING p.vec[6] = parms$ResVar.bolt[6] #Residula variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[6] #Average seedling size p.vec[13] = parms$size.dist.resse[6] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[6],parms$growth.intercept.nobolt.se[6]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[6],parms$growth.slope.nobolt.se[6]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[6],parms$growth.intercept.bolt.se[6]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[6] ,parms$growth.slope.bolt.se[6]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[6],parms$surv.intercept.se[6]) p.vec[8] = rnorm(1,parms$surv.slope[6],parms$surv.slope.se[6]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[6],parms$bolt.intercept.se[6]) p.vec[10] = rnorm(1,parms$bolt.slope[6],parms$bolt.slope.se[6]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[6],parms$recruit.se[6]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[6],parms$size.dist.unconditional.se[6]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[2],parms$zero.seed.intercept.se[2]) p.vec[15] = rnorm(1,parms$zero.seed.slope[2],parms$zero.seed.slope.se[2]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[2],parms$seed.intercept.se[2]) p.vec[17] = rnorm(1,parms$seed.slope[2],parms$seed.slope.se[2]) resBT6[n]=ParBoot(p.vec) } pick = (result$species =="BT") & (result$herb == "reduced")& (result$comp=="high") result[pick,]$lambda=median(resBT6) result[pick,]$LCI=quantile(resBT6,0.05) result[pick,]$UCI=quantile(resBT6,0.95) ########################################################################### ########################################################################### #TT, ambient herbivory, low competition resTT7=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[7] #Residula variance of regression line #GROWTH - BOLTING p.vec[6] = parms$ResVar.bolt[7] #Residula variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[7] #Average seedling size p.vec[13] = parms$size.dist.resse[7] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[7],parms$growth.intercept.nobolt.se[7]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[7],parms$growth.slope.nobolt.se[7]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[7],parms$growth.intercept.bolt.se[7]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[7] ,parms$growth.slope.bolt.se[7]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[7],parms$surv.intercept.se[7]) p.vec[8] = rnorm(1,parms$surv.slope[7],parms$surv.slope.se[7]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[7],parms$bolt.intercept.se[7]) p.vec[10] = rnorm(1,parms$bolt.slope[7],parms$bolt.slope.se[7]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[7],parms$recruit.se[7]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[7],parms$size.dist.unconditional.se[7]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[9],parms$zero.seed.intercept.se[9]) p.vec[15] = rnorm(1,parms$zero.seed.slope[9],parms$zero.seed.slope.se[9]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[9],parms$seed.intercept.se[9]) p.vec[17] = rnorm(1,parms$seed.slope[9],parms$seed.slope.se[9]) resTT7[n]=ParBoot(p.vec) } pick = (result$species =="TT") & (result$herb == "ambient")& (result$comp=="low") result[pick,]$lambda=median(resTT7) result[pick,]$LCI=quantile(resTT7,0.05) result[pick,]$UCI=quantile(resTT7,0.95) ########################################################## #TT, ambient herbivory, med competition resTT8=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[8] #Residula variance of regression line #GROWTH - BOLTING p.vec[6] = parms$ResVar.bolt[8] #Residula variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[8] #Average seedling size p.vec[13] = parms$size.dist.resse[8] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[8],parms$growth.intercept.nobolt.se[8]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[8],parms$growth.slope.nobolt.se[8]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[8],parms$growth.intercept.bolt.se[8]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[8] ,parms$growth.slope.bolt.se[8]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[8],parms$surv.intercept.se[8]) p.vec[8] = rnorm(1,parms$surv.slope[8],parms$surv.slope.se[8]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[8],parms$bolt.intercept.se[8]) p.vec[10] = rnorm(1,parms$bolt.slope[8],parms$bolt.slope.se[8]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[8],parms$recruit.se[8]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[8],parms$size.dist.unconditional.se[8]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[9],parms$zero.seed.intercept.se[9]) p.vec[15] = rnorm(1,parms$zero.seed.slope[9],parms$zero.seed.slope.se[9]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[9],parms$seed.intercept.se[9]) p.vec[17] = rnorm(1,parms$seed.slope[9],parms$seed.slope.se[9]) resTT8[n]=ParBoot(p.vec) } pick = (result$species =="TT") & (result$herb == "ambient")& (result$comp=="medium") result[pick,]$lambda=median(resTT8) result[pick,]$LCI=quantile(resTT8,0.05) result[pick,]$UCI=quantile(resTT8,0.95) ########################################################## #TT, ambient herbivory, high competition resTT9=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[9] #Residula variance of regression line #GROWTH - BOLTING p.vec[6] = parms$ResVar.bolt[9] #Residula variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[9] #Average seedling size p.vec[13] = parms$size.dist.resse[9] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[9],parms$growth.intercept.nobolt.se[9]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[9],parms$growth.slope.nobolt.se[9]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[9],parms$growth.intercept.bolt.se[9]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[9] ,parms$growth.slope.bolt.se[9]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[9],parms$surv.intercept.se[9]) p.vec[8] = rnorm(1,parms$surv.slope[9],parms$surv.slope.se[9]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[9],parms$bolt.intercept.se[9]) p.vec[10] = rnorm(1,parms$bolt.slope[9],parms$bolt.slope.se[9]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[9],parms$recruit.se[9]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[9],parms$size.dist.unconditional.se[9]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[9],parms$zero.seed.intercept.se[9]) p.vec[15] = rnorm(1,parms$zero.seed.slope[9],parms$zero.seed.slope.se[9]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[9],parms$seed.intercept.se[9]) p.vec[17] = rnorm(1,parms$seed.slope[9],parms$seed.slope.se[9]) resTT9[n]=ParBoot(p.vec) } pick = (result$species =="TT") & (result$herb == "ambient")& (result$comp=="high") result[pick,]$lambda=median(resTT9) result[pick,]$LCI=quantile(resTT9,0.05) result[pick,]$UCI=quantile(resTT9,0.95) ################################################################# #TT, reduced herbivory, low competition resTT10=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[10] #Residula variance of regression line #GROWTH - BOLTING p.vec[6] = parms$ResVar.bolt[10] #Residula variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[10] #Average seedling size p.vec[13] = parms$size.dist.resse[10] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[10],parms$growth.intercept.nobolt.se[10]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[10],parms$growth.slope.nobolt.se[10]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[10],parms$growth.intercept.bolt.se[10]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[10] ,parms$growth.slope.bolt.se[10]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[10],parms$surv.intercept.se[10]) p.vec[8] = rnorm(1,parms$surv.slope[10],parms$surv.slope.se[10]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[10],parms$bolt.intercept.se[10]) p.vec[10] = rnorm(1,parms$bolt.slope[10],parms$bolt.slope.se[10]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[10],parms$recruit.se[10]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[10],parms$size.dist.unconditional.se[10]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[12],parms$zero.seed.intercept.se[12]) p.vec[15] = rnorm(1,parms$zero.seed.slope[12],parms$zero.seed.slope.se[12]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[12],parms$seed.intercept.se[12]) p.vec[17] = rnorm(1,parms$seed.slope[12],parms$seed.slope.se[12]) resTT10[n]=ParBoot(p.vec) } pick = (result$species =="TT") & (result$herb == "reduced")& (result$comp=="low") result[pick,]$lambda=median(resTT10) result[pick,]$LCI=quantile(resTT10,0.05) result[pick,]$UCI=quantile(resTT10,0.95) #################################################################### #TT reduced herbivory medium competition resTT11=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[11] #Residula variance of regression line #GROWTH - BOLTING p.vec[6] = parms$ResVar.bolt[11] #Residula variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[11] #Average seedling size p.vec[13] = parms$size.dist.resse[11] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[11],parms$growth.intercept.nobolt.se[11]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[11],parms$growth.slope.nobolt.se[11]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[11],parms$growth.intercept.bolt.se[11]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[11] ,parms$growth.slope.bolt.se[11]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[11],parms$surv.intercept.se[11]) p.vec[8] = rnorm(1,parms$surv.slope[11],parms$surv.slope.se[11]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[11],parms$bolt.intercept.se[11]) p.vec[10] = rnorm(1,parms$bolt.slope[11],parms$bolt.slope.se[11]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[11],parms$recruit.se[11]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[11],parms$size.dist.unconditional.se[11]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[12],parms$zero.seed.intercept.se[12]) p.vec[15] = rnorm(1,parms$zero.seed.slope[12],parms$zero.seed.slope.se[12]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[12],parms$seed.intercept.se[12]) p.vec[17] = rnorm(1,parms$seed.slope[12],parms$seed.slope.se[12]) resTT11[n]=ParBoot(p.vec) } pick = (result$species =="TT") & (result$herb == "reduced")& (result$comp=="medium") result[pick,]$lambda=median(resTT11) result[pick,]$LCI=quantile(resTT11,0.05) result[pick,]$UCI=quantile(resTT11,0.95) ############################################################## #TT reduced herbivory, high competition resTT12=rep(0,Nrep) p.vec<-rep(0,17) #GROWTH - NO BOLTING p.vec[3] = parms$ResVar.nobolt[12] #Residula variance of regression line #GROWTH - BOLTING p.vec[12] = parms$ResVar.bolt[12] #Residula variance of regression line #SIZE DISTRIBUTION OF SEEDLINGS # biomass~N(N(avg,bb.ucse),bb.ucrse) #p.vec[12] = parms$size.dist.mean[12] #Average seedling size p.vec[13] = parms$size.dist.resse[12] #residual se (bb.ucrse) for (n in 1:Nrep) { #GROWTH - NO BOLTING p.vec[1] = rnorm(1,parms$growth.intercept.nobolt[12],parms$growth.intercept.nobolt.se[12]) #Growth Intercept p.vec[2] = rnorm(1,parms$growth.slope.nobolt[12],parms$growth.slope.nobolt.se[12]) #Growth Slope (no bolting) #GROWTH - BOLTING p.vec[4] = rnorm(1,parms$growth.intercept.bolt[12],parms$growth.intercept.bolt.se[12]) #Growth Intercept p.vec[5] = rnorm(1,parms$growth.slope.bolt[12] ,parms$growth.slope.bolt.se[12]) #Growth Slope (no bolting) #SURVIVAL p.vec[7] = rnorm(1,parms$surv.intercept[12],parms$surv.intercept.se[12]) p.vec[8] = rnorm(1,parms$surv.slope[12],parms$surv.slope.se[12]) #FLOWERING PROBABILITY p.vec[9] = rnorm(1,parms$bolt.intercept[12],parms$bolt.intercept.se[12]) p.vec[10] = rnorm(1,parms$bolt.slope[12],parms$bolt.slope.se[12]) #RECRUITMENT PROBABILITY (GERMINATING AND SURVIVING UNTIL THE NEXT POPULATION CENSUS) p.vec[11] = rnorm(1,parms$recruit.ave[12],parms$recruit.se[12]) #SIZE DISTRIBUTION OF SEEDLINGS p.vec[12] = rnorm(1, parms$size.dist.mean[12],parms$size.dist.unconditional.se[12]) #SEED PRODUCTION #PROBABILITY OF ZERO SEEDS p.vec[14] = rnorm(1,parms$zero.seed.intercept[2],parms$zero.seed.intercept.se[2]) p.vec[15] = rnorm(1,parms$zero.seed.slope[2],parms$zero.seed.slope.se[2]) #LOG SEEDS PRODUCED PER PLANT AS A FUNCTION OF PLANT SIZE p.vec[16] = rnorm(1,parms$seed.intercept[2],parms$seed.intercept.se[2]) p.vec[17] = rnorm(1,parms$seed.slope[2],parms$seed.slope.se[2]) resTT12[n]=ParBoot(p.vec) } pick = (result$species =="TT") & (result$herb == "reduced")& (result$comp=="high") result[pick,]$lambda=median(resTT12) result[pick,]$LCI=quantile(resTT12,0.05) result[pick,]$UCI=quantile(resTT12,0.95) library(gplots) barplot2(log(result[1:6,4]),names.arg=rep(c("Low","Med","High"),2), density=c(20,20,20,-1,-1,-1), ylab="",cex.lab=1.5, cex.axis=1.5,cex=1.5,ylim=c(-3,4.5), xlab="Degree of Interspecific Competition",main="Cirsium vulgare",cex.main=2,font.main = 4, plot.ci=TRUE,ci.l = log(result[1:6,]$LCI), ci.u = log(result[1:6,]$UCI)) mtext(expression(paste("log ", symbol(lambda))), side = 3,line=-1, cex=1.8, adj=-0.1) mtext("A", side = 3,line=2.5, cex=2, adj=-0.1) text(1.5,3.5,"Ambient H.",cex=1.5) text(6,3.5,"Reduced H.",cex=1.5) abline(h=0, lwd=2) barplot2(log(result[7:12,4]),names.arg=rep(c("Low","Med","High"),2), density=c(20,20,20,-1,-1,-1), ylab="",cex.lab=1.5, cex.axis=1.5,cex=1.5,ylim=c(-3,4.5), xlab="Degree of Interspecific Competition",main="Cirsium altissimum",cex.main=2,font.main = 4, plot.ci=TRUE,ci.l = log(result[7:12,]$LCI), ci.u = log(result[7:12,]$UCI)) mtext(expression(paste("log ", symbol(lambda))), side = 3,line=-1, cex=1.8, adj=-0.1) mtext("B", side = 3,line=2.5, cex=2, adj=-0.1) text(1.5,3.5,"Ambient H.",cex=1.5) text(6,3.5,"Reduced H.",cex=1.5) abline(h=0, lwd=2) ######################################################################## ######################################################################## ParBoot=function (p.vec){ # minimum (0.9*minimum size from data) and maximum sizes (1.1*maximum size from data) minsize=-3.5*0.8; maxsize=1.2*6.7; #============================================================================================# # (II) Demographic functions for computing the kernel #============================================================================================# # survival function s(x) sx<-function(x,params) { xbeta<-params[7]+x*params[8]; return(exp(xbeta)/(1+exp(xbeta))); } #flowering probability fx<- function(x,params) { xbeta<-params[9]+x*params[10]; return(exp(xbeta)/(1+exp(xbeta))); } # growth function g(x,y,params) no.bolt.gxy.f<-function(x,y,params) { mux<-params[1]+ x*params[2]; dnorm(y, mean=mux, sd=sqrt(params[3]), log = FALSE); } # growth function g(x, params) bolt.gxy.f<-function(x,y,params) { mux<-params[4]+ x*params[5]; dnorm(y, mean=mux, sd=sqrt(params[6]), log = FALSE); } #seedling production seeds<-function(x,params){ #calculate if any seeds are produced xbeta<-params[14]+x*params[15]; if (params[14]>998) {s=1 }else {s = exp(xbeta)/(1+exp(xbeta))}; # given you produce seeds, how many do you produce n=exp(params[16]+params[17]*x) # probability of germinating and surviving until the next population sensus = recruitment r=exp(params[11])/(1+exp(params[11])) #print(c(s, n, r)) return (s*n*r) } # combine s and g into p(x,y,params) for non-flowering plants pxy<-function(x,y,params) { g = no.bolt.gxy.f(x,y,params)*(1-fx(x,params))*sx(x,params) return (g) } # combine s and g into p(x,y,params) for flowering plants pfxy <- function(x,y,params){ g = bolt.gxy.f(x,y,params)*fx(x,params)*sx(x,params) return (g) } # fecundity function f(x,y,params) including number and size distribution of kids fxy<-function(x,y,params) { #plants need to survive, annd flower to produce seeds nkids <- seeds(x,params) #pfxy(x,y,params)* #include size distribution return(nkids*dnorm(y,mean=params[12],sd=params[13])); } #============================================================================================# # (III) Construct the component matrices and their transposes #============================================================================================# # Global variables for midpoint rule approximation # n.big.matrix is the number of mesh points for size, n.blow is the number of blowouts. n.big.matrix = 1000; ##### Put all component matrices into 3-dimensional arrays P<-array(NA,dim=c(n.big.matrix,n.big.matrix)) #P[j,i,a] will be h*P_{a-1}(x_j,x_i) B<-array(NA,dim=c(n.big.matrix,n.big.matrix)) #B[j,i,a] will be h*F_{a-1}(x_j,x_i) C<-array(NA,dim=c(n.big.matrix,n.big.matrix)) #B[j,i,a] will be h*F_{a-1}(x_j,x_i) L= minsize; U= maxsize; n = n.big.matrix # boundary points b and mesh points y b = L+c(0:n)*(U-L)/n; y = 0.5*(b[1:n]+b[2:(n+1)]); # step size for midpoint rule, see equations 4 and 5 h = y[2]-y[1] #Create P and B matrices # Note the transpose is used so we get P(y,x) and B(y,x) from functions that # compute p(x,y) and f(x,y). # Note the use of outer to compute each component matrix without looping. P<-h*t(outer(y,y,pxy,params=p.vec)) B<-h*t(outer(y,y,pfxy,params=p.vec)) C<-h*t(outer(y,y,fxy,params=p.vec)) # Create all the transposes (used when finding reproductive value and left e-vector) tP=P; tB=B; tP[,]=t(P[,]); tB[,]=t(B[,]); #tC = t(C) tC=C; tC[,]=t(C[,]) #============================================================================================# # (IV) Model iteration functions #============================================================================================# Nt=matrix(0,n.big.matrix); Nt1=Nt; # population now and next year #initializing population start.N = rep(1,n.big.matrix) Nt=start.N iteration=function(Nt) { Nt1=(P%*% Nt)+ C %*% (B%*% Nt); return(Nt1) } iteration.t=function(Nt) { Nt1=tB%*%(tC%*%Nt)+tP%*%Nt; return(Nt1) } #============================================================================================# # (V) Start using the model #============================================================================================# ##### Estimate lambda and w by iterating unperturbed matrix Nt=matrix(1,n.big.matrix); qmax=100; lam=1; tol=1.e-8; while(qmax>tol) { Nt1=iteration(Nt); qmax=sum(abs(Nt1-lam*Nt)); lam=sum(Nt1); Nt=Nt1/lam; # cat(lam,qmax,"\n"); } stable.dist=Nt/sum(Nt); lam.stable=lam; ##### Estimate lambda and v by transpose iteration Nt=matrix(1,n.big.matrix); qmax=1000; lam=1; tol=1.e-8; while(qmax>tol) { Nt1=iteration.t(Nt); qmax=sum(abs(Nt1-lam*Nt)); lam=sum(Nt1); Nt=Nt1/lam; # cat(lam,qmax,"\n"); } repro.val=Nt/sum(Nt); lam.stable.t=lam; return(lam.stable) }